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ABSTRACT 
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We perform a linear stability analysis of the axisymmetric, relativistic, self- 
q \ similar isothermal disk against non-axisymmetric perturbations. Two sets of 

neutral modes are discovered. The first set corresponds to marginally unstable 
C^i perturbations driven by gravitational radiation, and the other signals the onset of 

bifurcation to non-axisymmetric equilibrium solutions to the Einstein equations. 

£j 

> : 

■ 1. Introduction 

Cai & Shu (2002) put forward the hypothesis that relativistic disks may ultimately play 
as important a role in astrophysics as their spherical counterparts. There is evidence that 
supermassive blackholes in quasars were formed during the early stages of galaxy formation, 
where the major contributor to the mass-energy density was gas rather than stars. The 
dynamics of gas is highly dissipative, so it is easier for the gas to lose significant amounts 
of energy to reach the desired compactness for general relativistic effects to be important. 
More problematic is how to get rid of excess angular momentum if it is present initially. 

The pioneering work of Bardeen & Wagoner (1971) on uniformly rotating disks showed 
that with some assumed angular momentum loss, a Kerr blackhole may result in the collapse 
of such disks. However, as Mestel (1963) pointed out in a Newtonian context, there are 
astrophysical reasons to think that a disk specified by constant linear rotation velocity v is 
more realistic than one which has constant angular velocity. Through bar formation and 
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spiral density waves (for a review, e.g., see Shu et. al. 2000), such differentially rotating 
disks possess natural mechanisms for the outward transport of angular momentum and the 
inward transport of mass that would promote, in the relativistic regime, blackhole formation 
at the center. 

The first step in a systematic theoretical study of this possibility is to construct fully rel- 
ativistic, self-similar, rotating, flattened solutions. Lynden-Bell & Pineault (1978) performed 
such a pioneering analysis, but only in the cold limit. To suppress well-known axisymmetric 
instabilities that would fragment the disk into rings, it is necessary to include partial sup- 
port from isothermal pressure. If the pressure is exerted isotropically in three directions, the 
result (Cai & Shu 2003) is the relativistic generalization of the singular isothermal toroids 
(SITs) found by Toomre (1982) and Hayashi et. al. (1982). A useful approximation when 
such states become sufficiently flattened by rotation is to ignore the thermal dispersive speed 
in the vertical direction, while retaining it in the horizontal directions. In this approxima- 
tion SITs become completely flattened, singular isothermal disks (SIDs), whose equilibrium 
properties in the relativistic regime were studied by Cai & Shu (2002). These solutions are 
infinite in extent, possess infinite total mass, and contain a naked singularity at the origin 
(a "baby blackhole" with vanishingly small mass). 

The formal approximation of highly flattened configurations is satisfied for SITs only 
when the Mach number M = ^> 1. Nevertheless, , even if M ~ 1, Cai & Shu (2003) 

found that the critical condition under which the sequence of equilibria terminates as a 
function of M is nearly identical whether the pressure is exerted in three dimensions (SITs) 
or two (SIDs). The insensitivity of crucial properties of the equilibria to the assumption of 
infinitesimal thickness will hopefully carry over to the analysis of their stability. 

On dimensional grounds, if the disk becomes gravitationally unstable to overall gravita- 
tional collapse (a possibility if the disk is sufficiently slowly rotating), the mass of the baby 
black hole will grow linearly in time as a result of axisymmetric collapse. In the analogous 
problem of the collapse of relativistic singular isothermal sphere (SIS), Cai & Shu (2004) have 
shown that the growth of a black hole with finite mass introduces a (spherically symmetric) 
horizon which covers up the singularity. It is intriguing to ask whether such singularity 
in the case of a collapsing, relativistic SID will remain naked if the requirement of axial 
symmetry is relaxed. In order to answer this question, one must first construct relativistic 
SIDs that are non-axisymmetric states of equilibria. One of the goals of the present paper 
is to make a start on this problem, by finding the points of non-axisymmetric bifurcation 
along the sequence of axisymmetric SIDs, thereby generalizing the Newtonian work of Syer 
& Tremaine (1996) and Galli et. al. (2001). 

A feature with no Newtonian analog appears with the completion of the analysis: the 
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appearance of a secular instability which afflicts all relativistically rotating disks. This insta- 
bility, associated classically with the name Rossby (so the corresponding perturbations are 
called R- modes), arises because general relativity admits the radiation of angular momentum 
(and energy) by gravitational waves (Chandrasekhar 1970a,b). Basically, if the rotation of 
the underlying axisymmetric state is high enough, a counter-rotating disturbance appears 
corotating to an inertial observer. Such disturbances have negative angular momentum den- 
sity in the local rest frame of the disk. Since gravitational radiation carries away positive 
angular momentum, the amplitude of the R-mode perturbation grows in time. 

The phenomenon renders, in some sense, all astrophysically rotating systems potentially 
unstable to spin-down on a gravitational-radiation time scale. Whether the R-mode insta- 
bility competes with the gravitational torques associated with barlike deformations or spiral 
density waves remains a problem for future study. The self-similar models and techniques 
used in the present paper are capable only of determining the criterion for the onset of secular 
instabilities, and not their growth rates and evolution into the nonlinear regime. Thus, mod- 
ifications are still required for application to astrophysically realistic circumstances, where 
the origin does not contain a singularity from the start, and where spacetime at infinity is 
flat. 

In this paper, we restrict our study of the stability of relativistic SIDs to non-axisymmetric 
perturbations with the same scale-free character as the equilibrium state (i.e., with the same 
power-law radial dependence). In the nomenclature of Syer & Tremaine (1996), we consider 
only aligned perturbations, and no spiral disturbances. In section 2, we review the basic 
properties of axisymmetric SIDs. In section 3, we develop the mathematical formulation of 
the stability analysis, including the metric and matter perturbation. In section 4, the equa- 
tions are solved in the Newtonian limit, and the result is compared to Shu et. al. (2000). In 
section 6, the perturbation equations are solved in the full relativistic context, and we offer 
physical interpretation of the results. 



2. Review of Axisymmetric Disk Solution 

Start out with a self-similar axisymmetric metric 

ds i = -r 2n e N dt 2 + r 2 e 2p - N {d<\> - r n ~ l e N ~ p Qdt) 2 + e z ~ N {dr 2 + r 2 d9 2 ), (2-1) 

where N, P, Q, Z are functions of 9 and n is a constant measuring the strength of the 
gravitational field. For numerical convenience, we have chosen the equatorial plane to be 
at some polar angle 9q, which is determined as an eigenvalue of the problem. The locally 
nonrotating observer (LNRO) defines an orthonormal tetrad frame analogous to the inertial 
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frame: 



3 (i) 



(r- n e- 1 * N ,r- 1 Qe 1 2 N - p ,0,0) 



11 = (0,r _1 e 



-1 ijV-P 



0,0) 



e (2 ;=(0,0,e^),0) 
e (3) '*=(0,0,0,r- 1 e ^- z )). 



(2-2) 



We look for solutions to the Einstein field equations with a disk matter source described by 
a constant linear rotation velocity and a two-dimensional isotropic pressure. In the frame of 
a LNRO, the stress-energy tensor is taken to be 

J(o)(0) = — — , -Iff 



(i)(i) 



1 -v 2 
+ ev 



■ (o)(i) 



1 -v 2 



(2-3) 



1 -w 2 

where £ oc 8(9 — 9 ) and p^ = p T = js. Define 



, T(2)(2) = Pr, 



e = (l + n)fl, ' =5§> 



1 + n 



A = 5(0 -e ). 

After some algebra, part of the Einstein equations are cast into a set of dynamic equations 



N" = -N'P' 



2n 
1 +n 



+ Q 2 F 2 + Q 2 



1 — n 
1 + n 



+ e 



1 + 27 + 
1-v 2 



A, 



P" = -P' 2 - 1 + e^A, 



Q" = -Q'P' + Q 



:i-q 2 ) 



1 — n 
1 + n 



+ (jv' - p') 2 - Q 2 ^ 2 



(2-4) 



e ( X +7) ^ ^ A ' 



z" = -z'p' + 2g 2 



1-v 2 

1 — n 
1 + n 



An 2 



+ Q 2 F 2 + 2s 



1 + n) 2 



f +7 
1 -v 2 



where F = N' + logQ' — -P'- The rest of the field equations and the equation of motion form 
a set of constraint equations 



Q(e )u(l - n)(l + 7) + 7 + v 2 - n(l + 7) = 0, 

Z > - -^-N' - Q 2 ^^F = 0, 
1+n 1+n 

4n 2 



QHF 2 



n 



1 + n 



(2-5) 



2 cot ez' - AT' 2 + 



(l+n)< 



= 0. 
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Due to the contracted Bianchi identity, two of these equations are redundant, which were 
used for a numerical consistency check. For detailed discussion on the properties of these 
solutions, the readers are invited to refer to Cai & Shu (2002). 



3. Perturbed Configuration 

We will use the Eulerian description for the non-axisymmetric modes. The perturbation 
in the metric is 



^9 [iv 



< 1. 



9aj3 

We define the change in the contravariant components of the metric as 
so that 1 

(<r + 6<r)(g vp + 6g vp ) = &> p + 0(h 2 ). 

Notice that h^ v is not a tensor with respect to the unperturbed metric, and its projection 
onto LNRO is not a Lorentz scalar. Thus the directional derivatives of /ir )(&) wm involve 
more than the usual Ricci rotation coefficients, which destroys the simplicity of the tetrad 
formalism. As a result, we will compute in the coordinate frame whenever derivatives are 
involved. However, the tetrad frame defined by (2-2) does offer a clean separation of r from 
the other coordinates, so we shall project our results onto LNRO after the derivatives have 
been taken. 

The change in the Ricci tensor reads (see, e.g., Wald 1984) 

8R»„ = - {h a ^. ua + h a u . m - h^ u . a a - h ;flu ), (3-1) 

where the raising and lowering of indices are done with the unperturbed metric. Instead of 
computing the Einstein tensor, we work out the trace-reversed stress-energy tensor. Taking 
the direct variation of the stress-energy tensor, we have 

5(T, U - l -g^ v T) = 5T^ U - Ui^T - l -g^5T^g aP - l -g^K p . 

There is one subtlety in writing down this expression, due to the non-tensorial nature of the 
metric variations h^ u . Explicitly, 

8T^ = 5T^g a ,g Pv + T^h av + T v a h aiM ^ 5T^g aii g Pv . 



: It is unfortunate that we have two S symbols here - one denoting Eulerian change and the other denoting 
the Kronccker-Delta function. There should be little confusion from context, however. 
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To avoid such confusion, we shall adopt the convention that only contravariant components 
T Q/3 are varied 2 . So the Einstein field equation reads 

h a + h a — h a — h 

=8n(25T^g^ gfSl/ + 2T^h au + 2T v a h aix - h, u T - g^JT^g^ - g^T^h^). 



3.1. Gauge Choice 

Let's consider perturbations with time and angular dependence e Mn ^ >_Wi '*. On dimensional 
grounds, a scale-free disk can not support modes with uj ^ 0. The limiting case u = signals 
the onset of bifurcation or marginal stability of a particular mode. The most general form 
of ^(a)(6) ma Y be written as 



/1(a)(6) = 



^00 


h i 


h 2 


h 3^ 


hoi 


h u 


hi2 


hi 3 


h 2 


hi2 


h 2 2 


h 2 3 


\ho3 


hi3 


h 2 3 


h33/ 



where the 10 h entries are functions of 6 only. Geometrically, the metric coefficients are the 
inner products of the basis vectors 

d d 
^ dx^ dx u 

Since we are only considering polar perturbation, the system is symmetric about the equator 
and the metric is invariant under the diffeomorphism 9 — > 20 o — 6. This implies the boundary 
condition 

h t e = h^e = h r9 = =>- h^ = for // ^ 3 (3-3) 

on the disk. 

We proceed as follows. Project the left-hand side of (3-2) onto the LNRO, and write 
the result as Zr a )(6) (1 + n) 2 e irn ^ e N ~ z / r 2 . Expand l( a )(b) an d replace the second derivatives 
of zeroth-order metric coefficients with the unperturbed Einstein equations (2-4). This will 
introduce singular terms on the disk. Since we require the metric to be continuous across 
the disk, all singular terms must balance for the first-order equations in l< a )(b), which are L a \ 3 
(the second-order equations are acceptable since the first derivatives will in general have a 
jump there). Miraculously, with the condition (3-3), all first-order equations are regular. 



2 This is not entirely unfamiliar. Recall that in the super-Hamiltonian formalism, the conjugate momentum 
to x M is p^, which is what we vary, not pP. 
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To proceed further, we need to choose a gauge. Consider an infinitesimal coordinate 
transformation x M — > = + £ M (a;), where £ M is of the same magnitude as h^ v . This 
induces a transformation on the metric in the usual way, 



Thus, the coordinate freedom we have in general relativity corresponds to the gauge freedom 
h^ v + 2£( /1 .„). As suggested by the boundary condition on the disk, we shall promote 

(3-3) to a gauge condition. There is one more degree of freedom which we will fix here. The 
total gauge thus reads 

^03 = h 13 = h 23 = 0, /in = ^33- (3-4) 

The last condition resembles the Regge- Wheeler gauge in spherical symmetry. With the 
gauge condition, we may write 



h tt = r 2n e N [a + Q 2 b - 2Qd]e im<t> , 
h tr =tr n e z ^(f-Qj)e m t h 4 



h 



'be' 



t(j> = r-^e (d - Qb)e 



r 2 e 2P-N u Jm4 



im<j> 

i 

2P-Z h 

nee, 



V = ire p+z / 2 - N je im4 



(3-5) 



which corresponds to 



h 



(a)(6) 



= e 



The left-hand side of (3-2) now reads 



fa d if 0\ 

d b ij 

if ij c 

\0 bj 



(3-6) 



l 00 = - a " - (ijV' + P') a ' + \n'c' + 2QFd' 




eA \b 
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z/2 _ P 2Qm(l + 2n) z/2 „ p 2m(n + nQ 2 - Q 2 ) . 

(1 + n) 2 ; + 6 (1 + n) 2 h 

l u = -b" + (P' - \n')cl' - P'b' + ^-N' - P')c' + 2QFd' 

Z Z 



+ 



Z-2P 



m 



(1 + n)- 



Q 2 



1 — n 
1 + n 



+ F 2 



+ <e 



z _ 2P (l-Q 2 )m 2 ^ 2 



(1 + n) 



- Q 



+ <e 



Z-2P 



m 



(l+n)> 



+ g 2 



+ 2Q{ e 



Z-2P 



m 



(1 + n)' 



1 — n 
1 + n 

1 — n 
1 + n 



1 — n 
1 + n 

2 



F 



1 + n 



1 +n 



(N'-P')F \d 



+ 2Qe 



Z/2-P _ 



m 



-J - 2e z ' 2 - p - 



m 



-[Q 2 (l-n)+n + 2]j, 



(1 + n) 2 ' (1 + n 

I22 = -c" + \(Z'- N')a> + l -(N'-Z>- 2P')c' + Q 2 (j^f « 



+ 2Q 2 



1 — n\ 2 2n(l — n) 



1 + n 



+ Ul-Q 2 )e 



2\Z-2P 



m 



(1+n) 2 

2 



eA \b-2Q 



1 — n 



1 +n 



(1 + n)' 



+ eA-Q 2 



1 — n 
1 + n 



2n(l -n) 
(1 + n) 2 



+ Qe z/2-P_^_ f + eZ/2 , P _2m_ iQ2 _ 1 _ ng2)j . 



(1 + n)' 



hi = ~d" + \QFa' + QFb' + )-QFc' - P'd 

Z Z 

f 79P m 2 n 1 — n 
+ Q{ e z ~ 2P - - 2 



(1+n) 2 (1 + n) 2 

2 



b + Q{e 



,Z-2P 



m 



+ 2- 



1 — n 



(1+n) 2 (1 + n) 2 



+ {(i-Q 



1 — n 



-eA 1 + (N'-P') 2 -Q 2 F 2 }d 



- e z ' 2 - p 



l + n / 

2mn „ n _ 7/9 P m(n — 2) . 
f + 2Qe z ' 2 - p ■ J j, 



(1 + n)' 



(i + n y 



il 



02 



f" + P'f'-QFf+\Qe 



Z-2P 



m 



(i + n y 



QF + 2+t2tiA }3 



+ e 



Z/2-P _ 



m 



(1 + n)' 



1 1 

-Q(l - n)a + Q(l - n)b - -(n + 3)Qc - (1 - n)d 
z\ z\ 
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/ „Z-2P ™ 2 (M , 1 7 A2 , 2n 1 r) 2 F 2 | ^+7 ~ A 1 j 

« u = / - g*r + + g{-e*- ^ - 2^ + ( I* - N> )F } f 

+ ^""(T^ {(! - "X 1 + ^> + 2 (! - n ^ b 

- [1 + n + ^g 2 (l - n)]c - 3(1 - n)Qd| 

f 7 2P „2 m 2 ^f^~ n V 1 + 2n — n 2 
[ (1 + n) 2 \l + n/ (1+n) 2 

- (P> - Z'/2f + W - ^^~aU 

2 1 — v J 



3.2. Matter Content 

Recall that the disk is made of a two-dimensional perfect fluid. Explicitly, if we choose 
the equation of state p = je, the unperturbed stress-energy tensor may be written as 

= e[(l + 7 )u'V + 70""], for »,v = t, 0, r, and T Ae = 0. (3-7) 

In the presence of a perturbation, we still need to impose the condition that momentum flux 
and stress in the vertical direction vanish. Hence the first-order change in the stress-energy 
tensor is only for the upper-left 3x3 block: 

ST^ = Se[(l + 7 )«V + 70""] + e[(l + -f)(5u»u u + u»5u v ) - 7/1""]. 

As usual, the four-velocity is normalized, 

(«" + 5u»){u v + 5u v )( gia/ + V) = -!=*► = —u^h^. 

Projecting onto the LNRO frame, we have 

$ T (a)(b) = Se[(l + 7)u( fl )U(6) + 7"(a)(&)] + + 7)(<*«(a)U(6) + W(o)*«(6)) - 2 7^(a)(b)] (3-8) 
and 

5u^u (a) = -u^u^w (am , 

where 



1 „ ™„im0 .pinup o-yf 1 ' ■ 

= ( 1 _ u n n), faW = (_ p ye e ,o) 

\/l — f 2 \/l — f 2 Vl — v 2 \/l — v 2 \/l — V 2 
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With this parameterization, the normalization condition reads 

x = {y + d)v + Ua + v%). (3-9) 

Next, we work out the equation of motion (EOM) for the perturbed quantities. Although it 
is not needed if we solve all six unknown metric perturbations directly, the EOM provides a 
consistency check for the algebraic mess that is notorious in general relativity. Furthermore, 
as we will see later, for the self-similar disk, the EOMs are all algebraic, which are much 
easier to handle than the full Einstein equations. A direct variation of J 7 *" = reads 

5T» V . V + 5T» pv T pu + 5T u pu T^ = 0, 

where 

vp 2" \' L av;p ~ > L ap\v ,L pv;a)- 

Just as in the unperturbed EOM, the /i — 9 component needs to be satisfied identically. 
Physically, this component gives the fluid evolution in the 9 direction, which is trivial in this 
case. If we assume that the metric coefficients are even about the equatorial plane, then 
the first derivatives are odd, and vanish upon integrating through the plane. The resulting 
equation only contains h u g where v 7^ 3. This is another reason why the full gauge has 
to satisfy the boundary condition (3-3). Next, we consider the density perturbations. The 
zeroth-order density that is self-similar may be written as 

e = ±6(6-0 o ), 

where A = ee N °~ z °/S7r is some constant. Conventionally, the equatorial plane is located at 
9 = it/ 2 in spherical polar coordinates. However, to reduce eigenvalues of the problem, we 
rescaled 9 so that the disk is located at 6*o- When the geometry is perturbed, this "hidden" 
eigenvalue will in general change, hence we need to vary Q Q as well. Thus, we obtain 

X A im<f> A 

5s = —^5(9 - 9 ) + -M9 - O - 0id m *) - 6(9 - 9 )} 

(3-10) 

= e —[5A5{9 - 9 ) - AS' (6 - 9^]. 
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Projecting the right-hand side of (3-2) onto the LNRO and writing the result as if( a )( 6 )(l + 
n y e tm<f> e N-z j r 2 ^ -j.^ nonzero components are 3 



H 



oo 



1 — v z 



1 +7*; 2 + Qv(l +7) u(2u-Q)(l + 7) 



1 -w 2 



a + 



1 -v 2 



H 



01 



iH 02 
H n 

iH 12 
H22 

where 



1 -v 2 

1 — v z 
o (l + 7 )(l+^ 2 ) / 

^ 1- V 2 ^ 

7 + -fv 2 + 2 



1 -v 2 



£~A, 



(2 + t; 2 )(l +7) ^ 2(l + 7 t; 2 ) +7 + ^ 
u ~ _ tz 



1 -w 2 



1 -V 2 



£~A, 



1 - u 

,2 1 r» ,2 



1+7 1+7 



If 



(3-11) 



1 + v 2 + 27V' 



1 -w 2 



CA + 



7 + 3v 2 + 27W 5 
l-v 2 



b + 3v 



1 + 7 
l-v 2 



d + 4v 



1+7 
l-v 2 



£~A, 



-eA 



1 + 7 , . 7 + W + 2v 2 . I + 7 
2v- -/ H -„ j + 2v- -z 



l-v 

CA + 7 ceA, 



1 -v 2 



l-v 2 



CA = 87rfc 



r 2 e Z-N 

1 + n ' 



If we use the definition of 5e, the above expression simplifies to 

CA = £iA - ee Ar °- z °e z - JV A , ei = £~iA + e(Z - A^A =s> C = £i + e(Z' - N^Q,, 

where we integrated by parts on the last term. In fact, the explicit form of C is not required 
here since e.\ and Gi never appear independently in the equations. This observation suggests 
that Gi is a second order effect. 



3.3. Boundary Conditions 

/,From symmetry, all the non-axisymmetric metric components should vanish on the 
axis (where G = 0). Thus we can impose the conditions 



a = b = c = d = f = j = 0. 



(3-12) 



3 Actually #33 is nonzero too, but it contains exactly the singular terms from £33. 
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On the disk, the delta functions give rise to a jump in the derivatives of the perturbation 
functions using the second-order equations. Integrating across the disk, we have 

v 2 + 27 - ^v 2 - Qv(l + 7) 1 + 27 + 3v 2 + 2-fv 2 - Qv(l + 7) 



2d = e 

+ 
2b' = 



, -a + 

1 — v z 

(1 + -f)(5v - Q - Qv 2 ) 4u(l + 7 ) 



-d + 

1 — v z 1 — v z 

( 2^v 2 + 3v 2 + ^(1 + 7) 



-y 



1-v 2 

1 + 27 + v 2 



C, 



V 



1 > - b+ 1 

1 — v z 1 — V 



1-v 2 

, , . 1 + v 2 + 2-fv 2 



2c' = e{6-(l-7)c} + C, 
2d' = 



v(v 2 + 2)(l + 7 ), 37u 2 + 2 + u 2 , 2(l + v 2 )(l + 7 ) 
: s d : 7. V 



1-v 2 



1-v 2 



1-v 2 



2u(l + 7) 



c, 



= £ 



-2/' 
-2j" = e 



1-v 2 
2 + W -v< 



f 



2(1 



2u(l+7) l-7-2v 2 . 2u(l + 7) 



-/ + 



1 — v 2 1 — v 2 

The equations of motion in the disk read 



J - 



1-v 2 



2va 



v(v 2 + 2) + Q 



47V 2 + 7 + 2 + 3v 2 
1 + 7 



b-(v + Q)c - AQvd 



- 2(1 + v 2 + 2Qv)y e p - z/2 [2n - Qv(l - n)]z - 2 



m 



Q— — — + v 

1 + 7 



0, 



(1 +v 2 )a 



7(1 -v 2 ) 
1+7 



+ AQv + Qv 3 + 3v 2 }b-(v + Q)vc - 2(1 + v 2 )Qd 



2v 



- 2(Q + Qt; 2 + 2% e p ~ z/2 (l + n)z - 2^Qv + 



u 2 + 7 U 



1 +7 J £ 



vQa + (2Q + Qv 2 + 2v)vb + 2(Q + Qv 2 + u)d + 2m 



Q + v 



.Z/2-P 



/ 



+ 2(Q + Qt; 2 + 2% + 2wn^±V/ 2 - p j + 2m^e z / 2 - ? 2 = 0. 

1 — n 1 — n 

Equations (3-12), (3-13) and (3-14) form the complete set of boundary conditions. 



(3-13) 



(3-14) 



4. Newtonian Limit 



As shown in Shu et. al. (2000), the Newtonian magnetized isothermal singular disks 
allow bifurcation to non-axisymmetric equilibria if the rotation velocity is high enough. 
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Specifically, the onset of bifurcation occurs at 



777 

D 2 = , for m > 2, (4-1) 

where D is the ratio of rotation speed to magnetosonic speed (which equals the sound speed 
at zero magnetization). Our equation should yield the same result in the limit t/ < 1 and 
7 <sC 1. As usual, the Newtonian limit is recovered by setting 

g 00 = -l-2$ + 0(v 4 ), g 0j =O(v 3 ), gi^^ + Oiy 2 ), 

where $ is the Newtonian gravitational potential and is of order v 2 . We start with the 
unperturbed axisymmetric solution described by the metric (2-1) and the associated Einstein 
equations (2-4) and (2-5). The Newtonian limit becomes 



(4-2) 



or 



$ = (N/2 + nlogr)(l - Q 2 ) = 0{v 2 ) + 0(v 4 ), 
e 2p - N = sm 2 6 + 0(v 2 ), e p Q = 0(v 3 ), Z-N = 0(v 2 ), 

iV ~ n ~ Z ~ O ~e p — sin 9 ~ 7 ~ t> 2 , and Q ~ v 3 . 

The equation for iV reads 

N" + N' cot + 2n = eA, iV'(O) = 0, 

which has the solution 

iV' = — 2ntan— , 4ntan— = e. 

2 ' 2 

Since 6 ~ f , we have e 4n, which means it's also small. Thus, the equation for P 
becomes 

= 2 cot 6 =$> ©o = -, and = 6. 
Z may be most directly computed through the second constraint equation relating it to N', 

Z' = 2nN' = -An 2 tan °- = 0(v 4 ). 

Finally, the last constraint equation may be solved order by order. The 0(v 4 ) terms are 
identically zero by our solution of N'. The next order is 0(v 6 ), which reads 

(Q cot 9 - Q'f -Q 2 -2Q cot 6(Q cot 9 - Q r ) = 

=>• log Q' = esc 9 
q 

=^Q = Ctan -. 
^ 2 
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Integrating the dynamic equation for Q, we obtain a jump condition, which in this limit 

reads 

C = 2n{C + 2v) => C = Anv. 
Putting everything together, the limiting solution is 

n = v 2 + 7, e — 4n, P = log sin 9, 

^ A $ „ 2 ( 4 - 3 ) 

iV = -2ncot -, Q = 4ra;tan-, Z=-4rTtan-. 

2 2 2 

Of course, in the purely Newtonian case, Q and Z are taken to be since they are of higher 
order. A simple integration yields 

N = -2nlog(l - cos#)/2 $ = -(v 2 + 7) log ^(1 -cos0) 

- _ 

which is the correct result for a hot Mestel disk. 

In the presence of a perturbation, we still demand the Newtonian limit to be valid. 
Thus, the only nontrivial metric perturbation is in g o which is a, and 

g m = -r 2n e N (l-ae m *). 

Expanding everything to leading order in v, the Einstein equations reduce to 

2 

in 

-a" - cot do! + — = CA. (4-4) 
sin 2 6> 

This is nothing more than the Poisson equation for the perturbed potential V = —a/2. 

Next, we'll derive the Newtonian version of the equations of motion. From the Poisson 
equation, we know that a is of order (, which is in turn of order v 2 , while y and z are both 
of order v. It is worthwhile to point out that we are doing a two-parameter expansion, one 
in v and one in the perturbation. In this limit, the EOM becomes 

y + vi = 0, 

+ 7)4 = 0, <«> 
2 m e 

2y + mz = 0, 

where a is evaluated on the disk, of course. A not so trivial calculation by Galli et. al. 
(2001) shows that a = (/2m. Combining all three equations, we have (recall e = A(v 2 + 7)) 

v 2 + 7 2 2v 2 v 2 m 

h v —7 = 0=^ — = or m = 1 . 

m m 2 7 m + 2 



- 15 - 



This is the Newtonian bifurcation point obtained by Shu et. al. (2000). 4 



5. Numerical Implementation 

Equation (3-2) is linear in the metric perturbations. Therefore it is ideal for finite 
differencing, which transforms the differential equations into a matrix equation. Consider a 
vector 

V = a©b©c©d©f©j©x, (5-1) 

where 

a = (ai, a 2 , ...a N ), a { = a(Oj), = 6 max , 

etc., and 

x = (e 1 ,y,iz). 

We do not include the values of the perturbation on the pole, since they all vanish by 
the boundary conditions (3-12). For simplicity, let's use a uniform grid as we did in the 
unperturbed solution. Thus 

Qi =ih, h = ma x/iV. 
With these definitions, the differential operators may be written as 

tti = 2h ' Qi = V • (5 " 2) 

It is easy to check that these expressions are accurate to second order. On the boundary, 
the situation is a little bit trickier. To find the correct differencing scheme, let's expand the 
function on the axis to second order: 

ai = a + ha' + ^h 2 ciQ, 
d2 = a + 2hd + 2h 2 aQ. 

Taking the proper linear combinations, we have 

, 4ai - 3a - &2 

9 2/ \ (5-3) 
„ _ a 2 - 2ai + a 

%~ ^ • 



4 The conclusion that eccentric m = 1 bifurcations can occur at any level of disk rotation is flawed, as are 
the analyses of Syer & Tremaine (1996); Shu et. al. (2000) and Galli et. al. (2001) on this point (Toomre, 
personal communication; Shu, Toomre, Tremaine, in preparation). Except for one special rotation rate (the 
true bifurcation value), the stress tensor is nonzero at the origin (implying a physically unrealistic steady 
flow of momentum from the origin to infinity), even though it has zero divergence. 
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Similarly, on the disk, 



1 

-i 

2 

a N _ 2 = a N — 2ha' N + 2h 2 a" N . 



Q-N-i — o-n — ha' N + —h 2 a'l^, 



Thus, 

, 3ajv — 4ajv„i + ajv_2 

o i_ (5-4) 

// _ a N-2 — jO-N-l + ON 

The set of equations (3-2) may now be written as a matrix equation 

M • V = 0. (5-5) 

It is not too hard to convince oneself that the matrix M is (6iV + 3) x (6iV + 3). After 
finite differencing, each component of kj is represented by a (67V + 3) x N submatrix, where 
the left-hand side is evaluated at = h,2h, ..(N — l)h. The last row in this submatrix is 
replaced by the boundary condition (3-13). We thus fill the first 6iV rows of M. The very 
last three rows are the remaining boundary conditions of (3-12). 

In order for (5-5) to have nontrivial solutions V, the determinant of M must vanish. 
Schematically, the elements of M are nonlinear functions of the azimuthal quantum number 
m, and unperturbed metric coefficients (which are, in turn, functions of v and 7). Thus, for 
a given value of m and 7, the problem of stability analysis is reduced to root finding of the 
equation 

|M m »|=0. (5-6) 



6. Results 

We scan the entire solution space looking for the solution to (5-6). As expected, the 
bifurcation points form two sets of tracks in the 7 — v 2 space. The behaviors of these two 
tracks are fundamentally different, since they are caused by two different mechanisms. We 
will discuss them separately. 



6.1. Radiation Driven Neutral Modes 

The first set of tracks is shown in Fig 1. These modes are believed to be the analog of 
Rossby modes first discovered by Chandrasekhar in 1970 and subsequently studied exten- 
sively in the context of neutron stars. Even though our self-similar disk geometry is in some 
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sense infinitely different from the finite spherical geometry of neutron stars, the underlying 
mechanism for these neutral modes can still be understood if one is comfortable with the 
idea of gravitational radiation with infinite wavelength. To make a bad situation even worse, 
since our disks are formally infinite in size, one is never able to reach the radiation zone to 
study the gravitational wave. However, since the original argument of Chandrasekhar did 
not rely crucially on the asymptotic flatness of spacetime, the qualitative result still holds 
in this pathological case. Consider a non-axisymmetric disturbance in the disk which moves 
at a velocity V\ < v. As a result, the total angular momentum (or more appropriately, the 
specific angular momentum) decreases. These perturbations will in general radiate due to 
non-axisymmetry. If in the LNRO V\ < 0, then gravitational radiation carries away nega- 
tive angular momentum, which damps the amplitude of perturbation, and these modes are 
stable. On the other hand, if v\ > 0, gravitational radiation carries away positive angular 
momentum, and thus the amplitude of perturbation has to grow to make the total angular 
momentum more negative. 

For a given equation of state specified by 7, the radiation-driven neutral modes occur 
at lower Mach number for increasing m. This is expected from the analysis of Friedman & 
Schutz (1975, 1978a,b). In fact, the onset of instability for m — > 00 occurs at zero rotational 
velocity. However, for a realistic system, the strength of a particular unstable mode is 
intimately related to the magnitude of the imaginary part of its frequency. For the m — > 00 
mode, even though it is formally unstable at zero rotation, the characteristic growth time 
scale is infinite. This is due to the fact that multipole radiation is exceedingly weak for 
higher values of m. If we truncate the self-similar disk, the strongest unstable modes - that 
is, the modes with shortest growth time - are still the ones with small m. In the absence of 
viscosity, these modes will grow until the non-linear effects set in and limit the final rotation 
speed of the full disk. 

In addition to the R-mode tracks, the Q — 1 curve is also plotted in Fig 1. We would like 
to remind the readers that these tracks represent models where the characteristic frequency 
of a given mode becomes purely real, i.e., modes that are marginally stable. It is not too 
difficult to convince oneself that every mode in the background of an ergoregion is unstable. 
Using the simplistic picture of retrograde disturbances, we see that in the ergoregion, every 
mode has to propagate in the direction of the underlying disk as seen by an LNRO, and grav- 
itational radiation will drive them unstable. Furthermore, as Friedman (1978) demonstrated, 
a spacetime with an ergoregion is unstable even under scalar and vector perturbations. From 
Fig 1, the m = 2 track is entirely above the Q — 1 curve and so is part of the m = 3 track. 
In the presence of perturbing matter fields other than those in the disk itself, the ergoregion 
is likely to put a more stringent limit on the maximum rotation rate for a disk. 
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The question arises whether the growth of R-modes might be suppressed by viscous 
torques due to, e.g., a magneto-rotational instability (MRI) acting in the disk Balbus & 
Hawley (1991). Since the viscous effects must act to erase non-axial symmetries of length 
scale ~ r/m, the damping time associated with linear global R-modes must be the diffusion 
time scale to ~ r 2 /(m 2 z/), where v is the effective MRI kinematic viscosity. If relativistic 
disks are even weakly magnetized and electrically conducting, as their Newtonian counter- 
parts are believed to be, for m ~ 1, the time scale to may be only two orders of magnitude 
longer than the dynamical timescale r/v. R-mode spin-down for such disks may then be 
effectively suppressed by the MRI viscosity, but calculations of non-self-similar relativistic 
disks are needed to answer definitively whether the growth rate of R-modes (here zero) can 
overcome the viscous damping. 



6.2. Newtonian Bifurcation Track 

The second set is the generalization of Newtonian bifurcation computed by Shu et. al. 
(2000). In Fig 2, we plotted these extended "Newtonian" tracks for m = 2,3,4,5, and 
oo. The finite-m values are the numerical result from solving (5-6). When m becomes 
large, we may approximate it by a continuous variable and use it as a parameter in the 
asymptotic expansion. Assume all the perturbation amplitudes remain infinitesimal in the 
large-m limit, then the coefficients of each power of m need to vanish independently. This 
requirement translates to 

a + 2Q 2 b + Q 2 c - 2d = 0, Q(l + 2n)f = (n + nQ 2 - Q 2 )j, 
- a + (1 - Q 2 )b + c + 2Qd = 0, Qf= [Q 2 {1 -n)+n + 2}j, 
c=0, nQf = -(Q 2 -l-nQ 2 )j, 

Qb + Qc=0, nf = Q(n-2)j, (6-1) 
Iq(1 _ n )a + Q(l - n )b-±(n + 3)Qc - (1 - n)d = 0, Qj = f, Qf = Q 2 j, 

(1 -n)(l + \Q 2 )a + 2(1 - n)Q 2 b - [1 + n + ]-Q 2 {l - n)]c - 3(1 - n)Qd. 

Regardless of whether the disk is rotating, these (linear, homogeneous) equations only 
have trivial solutions. This is the familiar Cowling approximation where the metric perturba- 
tions vanish in the large-m limit. These Cowling modes may be understood in the following 
schematic way. Recall that in the Newtonian limit, we can invert Poisson's equation via 
a Green's function, and obtain an integral representation of the gravitational field. This 
procedure can also be done for the full Einstein equations in principle, although not analyt- 
ically. Mathematically, the integral representation of metric perturbations is effectively an 
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average of the matter perturbation times the Green's function. Therefore, in the m — > oo 
limit, the gravitational field is indifferent to the matter perturbation, since it averages to the 
axisymmetric equilibrium over any finite angular integration. This line of argument may be 
made mathematically rigorous with some more thought, but it is not necessary here. 

In the absence of a metric perturbation, the equations of motion now simplify to 



(l + v 2 + 2Qv)y + 



1 



Q——- + v 

1 + 7 



4 = 0, (Q + Qv 2 + 2v)y + 

e 

Eliminating the factor (/ey, we can combine these equations to give 

(1 + v 2 + 2Qv)(Qv + Qvy + v 2 + 7) = (Q + Qjv 2 + v + vj)(Q + Qv 2 + 2v). 

The last equation implicitly defines a surface of neutral modes Q(9q) = 7(7, v). Recall that 
the solution space of axisymmetric disks can also be viewed as a two-dimensional surface 
defined by Q(6q) = (7(7, v). Thus the intersection of these two surfaces gives the bifurcation 
curve in the 7 — v space. This curve is also plotted in Fig. 2. 

Near the origin, these bifurcations tracks recover the Newtonian result, where the bi- 
furcation point is located at 

v 2 m 

— = -z, ^,7^0. 

7 m + 2 

As 7 increases, the relativistic effects become important. Intriguingly, the point of bifurcation 
occurs at lower Mach number for a relativistic disk than in its Newtonian counterpart. The 
curve defined by Q = 1 on the disk is again plotted. As seen in Fig 2, these tracks are 
confined in the portion of solution space where the time-like Killing vector remains time- 
like. In fact, the only place where a bifurcating neutral mode is allowed in the ergoregion is 
for m = 00 and 7 = 1. Even then, the spacetime is only marginally ergo-like in the sense that 
Q — > l - on the disk. To understand this phenomenon, we examine the velocity perturbation 
corresponding to the bifurcating neutral modes. In the Newtonian limit, we can evaluate it 
analytically [see equation (4-5), or equation (21) of Shu et. al. (2000)] 

5v = ve im<t> = 5e 

V 4(t; 2 + 7) 

In general, this component of the eigenvector needs to be computed numerically. Actually, 
the exact form, or even the magnitude of this velocity perturbation, is not important. It 
suffices to know that Sv for this mode always has the opposite sign of Se. In the linear 
perturbation regime, where we can still apply the superposition principle, this observation 
has the following physical picture. The full disk solution has two components. One is the 
axisymmetric equilibrium rotating at v in the +0 direction with energy density given by e oc 



Qv + ——- 

1 + 7 
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e5{6 — 9 )/r 2 . The second component is the non-axisymmetric perturbation, which is a disk 
of infinitesimal energy density Se, and infinitesimal velocity field 5u^. Along the bifurcation 
track, this perturbation disk is always counter-rotating. The empirical result of Fig 2 leads us 
to postulate that the existence of counter-rotating non-axisymmetric disturbance is another 
necessary condition for bifurcation, at least for the disk geometry. A corollary is that the 
full non-axisymmetric disk spacetime can not have a stable ergoregion. This conjecture 
leads naturally to the confinement of bifurcation tracks in the portion of the solution space 
without ergoregions. For models with Q > 1 on the disk, the time-like Killing vector becomes 
space-like, and thus no counter-rotating trajectory is allowed. 

Can we elevate this conjecture to apply to a more generic relativistic non-axisymmetric 
equilibrium? The answer is a cautious yes. Without digressing too much into the mathe- 
matical structure of Riemannian geometry, we would like to offer the following plausibility 
argument. Suppose we are able to construct a fully nonlinear, non-axisymmetric stationary 
solution to the Einstein field equations. It can not have an event horizon, since black holes 
can not have "hair" (in this case, "hair" refers to mass multipole moments). Furthermore, 
in the absence of a space-like Killing vector d^, the non-vanishing component of g t j will 
generate a time-dependent quadrupole moment as seen by an inertial observer. As a result, 
gravitational radiation will continue to carry away angular momentum and energy until the 
system is either axisymmetric or static. In this aspect, the relativistic non-axisymmetric 
equilibria are analogous to the Dedekind ellipsoids, where the figure axes are static in an 
inertial frame, and the configuration is supported by pressure and internal motion. It can 
be shown (see, e.g., Chandrasekhar 1983) that a static metric can always be brought to the 
diagonal form after appropriate coordinate transformations. With a diagonal metric, the 
ergoregion defined by gu = coincides with the event horizon. Therefore, the absence of an 
event horizon means a non-axisymmetric static solution does not have an ergoregion. 



7. Summary and Discussion 

We performed a linear stability analysis of the relativistic self-similar disk against non- 
axisymmetric perturbations. For simplicity, we restricted the class of perturbation under 
consideration to be self-similar and polar. Mathematically, this means the scaling law is pre- 
served and the metric is symmetric about the midplane even in the presence of perturbation. 

As expected, the Newtonian bifurcations found by Shu et. al. (2000) and Galli et. 
al. (2001) have extensions into the fully relativistic regime. These tracks seem to exist 
only in models which do not admit an ergoregion. The corresponding velocity perturbation 
is strictly negative for any positive energy density increase. We thus hypothesize that in 
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addition to the existence of neutral modes, retrograde disturbance may also be a necessary 
condition for bifurcation to non-axisymmetric disk equilibria. This line of arguments leads 
us to speculate that the non-axisymmetric equilibrium solutions in general can not have 
ergoregions. We have no proof that the behavior is generic, whereas the bifurcation of 
rapidly rotating axisymmetric equilibria to non-axisymmetric forms probably is. 

In addition, we have discovered the onset of R-mode instability, which is driven by 
gravitational radiation. The marginal stability tracks follow the qualitative behavior first 
discussed by Friedman & Schutz (1978b), and the result here is probably also generic. For 
a self-similar disk, the entire m = 2 neutral-mode tracks and part of m = 3 occur in models 
with ergocones. We believe that in general for a fixed equation of state, the onset of instability 
occurs either at the ergoregion formation, or the tracks we computed, whichever have lower 
velocity. 

These studies are the first step to constructing fully non-axisymmetric relativistic equi- 
libria. For a given value of 7, if the axisymmetric state contracts quasi-statically by shedding 
angular momentum, the linear rotational velocity will increase. This evolution represents 
a vertical line in Figs 1 and 2. Eventually, the velocity will reach a value where a non- 
axisymmetric mode becomes unstable. If they undergo gravitational collapse, the secular 
instability will most likely survive over many dynamic time scales (which, in the purely self- 
similar case, is infinite). Therefore, the collapse will be fundamentally non-axisymmetric. 
For simplicity, we have only considered each Fourier component independently. In the linear 
perturbation regime, the effect of a general non-axisymmetric disturbance can always be de- 
composed into its Fourier components with each component decoupled from others. When 
the amplitudes become finite, our Fourier series will fail to converge, and a more detailed 
analysis is required. 

The current work serves as a spring board to one of the ultimate challenges in numer- 
ical relativity - the fully nonlinear numerical simulation of a non-axisymmetric collapse. 
Only then can we answer questions such as whether the central singularities of objects like 
relativistic SIDs remain naked. 

We thank the referee for valuable comments. This work has been supported by Academia 
Sinica through the Distinguished Postdoctoral Research Fellowship awarded to MJC, and by 
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Fig. 1. — Radiation driven neutral modes. For large enough values of m, the onset of 
instability is believed to occur at infinitesimal velocities. 
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Fig. 2. — The neutral mode curves which are connected to the Newtonian bifurcation. In 
the lower-left corner, where 7^0 and v — > 0, the slope is m/(m + 2), as expected from 
the Newtonian limit. The dashed line separates the solution space into ergoregion and non- 
ergoregion. 



